#ifndef STATSxx_MACHINE_LEARNING_GAUSSIAN_PROCESS_GAUSSIANPROCESS_HPP
#define STATSxx_MACHINE_LEARNING_GAUSSIAN_PROCESS_GAUSSIANPROCESS_HPP


// STL
#include <memory>  // std::unique_ptr<>
#include <utility> // std::pair<>
#include <vector>  // std::vector<>

// Boost
#include <boost/serialization/export.hpp>
#include <boost/serialization/serialization.hpp>     // boost::serialization::
#include <boost/serialization/unique_ptr.hpp>        // serialize std::unique_ptr<>

// jScience
#include "jScience/linalg/Matrix.hpp"                // Matrix<>
#include "jScience/linalg/Vector.hpp"                // Vector<>

// stats++
#include "statsxx/statistics/CovarianceFunction.hpp" // CovarianceFunction
#include "statsxx/statistics/covariance_functions.hpp"   // squared exponential


namespace Gaussian_process
{
    //=========================================================
    // GAUSSIAN PROCESS
    //=========================================================
    class GaussianProcess
    {

    public:

        // NOTE: A Gaussian process is fully specified by a mean function mu(x) and covariance function k(x,x').
        GaussianProcess();
        GaussianProcess(
                        std::unique_ptr<CovarianceFunction> const &k_arg
                        );
        ~GaussianProcess();

        // adding data:
        void add_data(const Matrix<double> &X_arg,
                      const Vector<double> &y_arg,
                      const double sigma_n_arg);

        double likelihood();
        std::vector<double> dlikelihood();

        // prediction:
        double predict(const Vector<double> &x) const;
        Vector<double> predict_deriv(const Vector<double> &x);

        // parameters:
        // NOTE: This needs to be public for the _GP package.
        std::unique_ptr<CovarianceFunction> k;

    private:

        // training information:
        //    bool data_addded;

        Matrix<double> X;
        Vector<double> y;
        double         sigma_n;

        Matrix<double> Ky;    // Ky = K + sigma_n^2*I
        Matrix<double> Kyinv; // Ky^-1
        Vector<double> alpha; // alpha = Ky^-1*y

        /*
         double L_LOO();
         std::vector<double> dL_LOO();
         */

        // ----- (Boost) SERIALIZATION -----
        friend class boost::serialization::access;

        template<class Archive>
        void serialize(
                       Archive            &ar,
                       const unsigned int  version
                       )
        {
            //ar.template register_type<SquaredExponential>();

            //ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(*this->k);
            //ar & boost::serialization::base_object<CovarianceFunction>(*this->k);
            ar & this->k;

            ar & this->X;
            ar & this->y;
            ar & this->sigma_n;

            ar & this->Ky;
            ar & this->Kyinv;
            ar & this->alpha;
        }

    };
}

#include "statsxx/machine_learning/Gaussian_process/src/GaussianProcess.cpp"

//BOOST_CLASS_EXPORT_KEY(Gaussian_process::GaussianProcess);


#endif
